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Abstract 

We present a novel method for identifying a biochemical reaction network based on 
multiple sets of estimated reaction rates in the corresponding reaction rate equations 
arriving from various (possibly different) experiments. The current method, unlike 
some of the graphical approaches proposed in the literature, uses the values of the 
experimental measurements only relative to the geometry of the biochemical reactions 
under the assumption that the underlying reaction network is the same for all the 
experiments. The proposed approach utilizes algebraic statistical methods in order 
to parametrize the set of possible reactions so as to identify the most likely network 
structure, and is easily scalable to very complicated biochemical systems involving a 
large number of species and reactions. The method is illustrated with a numerical 
example of a hypothetical network arising form a "mass transfer" -type model. 

Keywords: Biochemical reaction network, law of mass action, algebraic statistical 
model, polyhedral geometry. 
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1 Introduction 

In modern biological research, it is very common to collect detailed information on time- 
dependent chemical concentration data for large networks of biochemical reactions (see 
survey papers j2[[TT]). Often, the main purpose of collecting such data is to identify 
the exact structure of a network of chemical reactions for which the identity of the 
chemical species present in the network is known but a priori no information is available 
on the species interactions. The problem is of interest both in the setting of classical 
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theoretical chemistry, as well as, more recently, in the context of molecular and systems 
biology problems and as such has received a lot of attention in the literature over last 
several decades as evidenced by multiple papers devoted to the topic [IJ EH III El [10l fl5l 

nana nig. 

In general, two very different reaction networks might generate identical mass-action 
dynamical system models, making it impossible to discriminate between them, even if 
one is given experimental data of perfect accuracy and unlimited temporal resolution. 
Sometimes this lack of uniqueness is referred to as the "fundamental dogma of chemical 
kinetics" , although it is actually not a well known fact in the biochemistry or chemical 
engineering communities [3J |3J H] . Necessary and sufficient conditions for two reaction 
networks to give rise to the same deterministic dynamical system model (i.e., the same 
reaction rate equations) are described in [I], where the problem of identifiability of 
reaction networks given high accuracy data was analyzed in detail. The key observation 
is that, if we think of reactions as vectors, it is possible for different sets of such vectors 
to span the same positive cones, or at least to span positive cones that have nonempty 
intersection (see Figure [l] for an example) . 

On the other hand, it is often the case that experimental measurements for the 
study of a specific reaction network or pathway are being collected under many different 
experimental conditions, which affect the values of reaction rate parameters. Almost 
always, the reactions of interest are not "elementary reactions" , for which the reaction 
rates parameters must be constant, but they are so called "overall reactions" that 
summarize several elementary reaction steps. In that case the reaction rates parameters 
may reflect the concentrations of biochemical species which have not been included 
explicitly in the model. In such circumstances the reaction rate parameters are not 
constant, but rather depend on specific experimental conditions, such as concentrations 
of enzymes and other intermediate species. Therefore, the estimated vector of reaction 
rate parameters will not be the same for all experimental conditions, but each specific 
experimental setting will give rise to one such vector of parameters. However, the set 
of all these vectors should span a specific cone, whose extreme rays should identify 
exactly the set of reactions that gave rise to the data. 

The purpose of the current paper is to propose a statistical method based on the 
above geometric considerations, which allows one to take advantage of the inherent 
stochasticity in the data, in order to determine the unique reaction network that can 
best account for the results of all the available experiments pooled together. The idea 
is related to the notion of an algebraic statistical model (as described in |12j Chapter 1), 
and relies on mapping the estimated reaction parameters into an appropriate convex 
region of the span of reaction vectors of a network, using the underlying geometry to 
identify the reactions which are most likely to span that region. As shown below, this 
approach reduces the network identification problem to a statistical inference problem 
for the parameters of a multinomial distribution, which may then be solved for instance 
using the classical likelihood methods. 
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Figure 1: Identifiability of reaction networks given experimental data. Note that, for a determinis- 
tic mass-action model, even if we can estimate the vector K of parameter values with great accuracy, 
we cannot determine if the "correct" reaction network is {Aq — > 2A\ : Aq — > A\ + A2, Aq — > 2A3} 
or {Aq — » 2A2,Ao —> A% + A^,Aq — » 2A3}, because if belongs to the span of either one of these 
networks. However, if instead a single point K one has available a set I? = {Ki,i = 1 . . . , k}, 
interpreted as a result of random selection of parameter rate values according to some probability 
law, then the spanning cone of the data points may be used to identify the sets of reactions that 
"best explain" the data. 

2 Maximum Likelihood Inference for a Biochemical 
Reaction Network 

In this section we develop a formal way of inferring a most likely subnetwork of a given 
conic network (i.e., network represented by a cone like the one in Figure [TJ of the 
minimal spanning dimension. For the inference purpose, in the network of m reactions 
we assume that the empirical data T> — {Ki,i — 1 . . . , k} c R d is available in the 
form of (multiple) estimates of the parameters of the system of differential equations 
corresponding to a hypothesized biochemical network. As illustrated in Figure [TJ such 
networks are in general "unidentifiable" in the sense that different chemical reaction 
networks may give rise to the same system of differential equations. However, in 
the stochastic or "statistical" sense it is possible to identify the "most likely" (i.e., 
maximizing the appropriate likelihood function) network as indicated by the data T>. 
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2.1 Multinomial model 



Consider d species, and m possible reactions with reaction vectors R = {R\, . . . , R m } C 
R d among the species. (For more details about how each reaction generates a reaction 
vector see pQ.) Let IZd denote the collection of all (™) positive cones spanned by 
subsets of d reactions in R. 

Denote by cone(R) the positive cone generated by the reaction vectors in R. Let 
S be the partition of cone(R) obtained by all possible intersections of non-dcgcncratc 
cones in TZd- Suppose S contains n full-dimensional regions S\, . . . , S n ; throughout we 
shall refer to these regions as building blocks, and to n as the number of building blocks. 

Let A m _! be a probability simplex in M. m and let 9 S A m _ 1 be a vector of prob- 
abilities associated with the reactions that give rise to R. We assume that these m 
reactions have the same source complex (i.e., form a conic network), since, as explained 
in [1], the identifiability of a network can be addressed one source complex at a time. 
Define the polynomial map 

g : A m _! -> R" 

where 



C=cone(R CT ( 1 ),...,R tr ( ti ))e7?. < j 

for i — 1 . . . , n. 

tioi(cns,) _ 



We take 



'^f- = if vol(C) = 0. Define = £ ff • • • 8 a(d) and 
p(^) = (pi(0) . . . ,p„(0)) = (gi(9)M6), . . .,g n (6)/s(6)). (2) 



In this setting p £ M. n is our statistical model for the data, after we substitute 
9 m = 1 — Sjli 1 Qj- Note that we may interpret the monomials 9 a (i) ' ' ' ^a(d) m @ 
as the probabilities of a given data point being generated by the d-tuple of reactions 
cr(l), . . . , a(d). With this interpretation the coordinate pt of the map p in ^ is sim- 
ply the conditional probability that the data point is observed in Si given that it was 
generated by a d-tuple of reactions. Note that the map p is rational but, as we shall 
see below, the model may be re-parametrized into an equivalent one involving only the 
multilinear map ( 1 ) . 

Let Ui denote the number of data points in Si . The log- likelihood function corre- 
sponding to a given data allocation is 

n 

l(6) = J2uilog P i(e). (3) 
1=1 

Our inference problem is to find 

rn 

6 = argmax e Z(0) subject to Oj = 1 and 9i > 0. (4) 



1 In general, it may be beneficial to consider various measures vol(-) which are absolutely continuous 
w.r.t. the usual Lebesque measure. For instance in Section 3 we describe an example where this measure 
is defined via gamma densities. 
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Example 2.1. Consider the two reaction networks described in Figure [7] Since the 
species Aq does not appear as a product of any reaction, the model has effectively 
d = 3 species and a total of m = 5 possible reactions R — {R\, ■ . ■ , R5} = {Aq — > 
2Ai,Aq — > Ai + A 2 , A — * 2A 3 , A — * 2A 2 ,A — ► Ai + A 3 }. In iftis case t/iere are 
n = 5 building blocks Si,...,Ss defined by the intersections of all non-trivial reaction 
cones generated by reaction triples. Thus denoting Cjkl = cone(Rj,Rk,Ri) for any 
triple {j, k, 1} € {1, . . . , 5} we have 

Si =Ci34 R C234 H C345 

5*2 =Ci34 n C145 n C234 n C245 

53 =Ci23 n C134 n C235 n C 3 45 

54 =Ci2 3 n C134 n C145 n C235 n C245 
5*5 =Ci23 n C125 n C134 n C145. 



Note that the cones C124 and C135 are degenerate and are not involved in the definitions 
of the Si 's. Denoting further v^ k \ = vol(Cjki fl Si)/vol(Cjki) for any triple {j, k, 1} € 
{1, . . . ,5}, we see that the the map (JlJ becomes 

gi {6) =v { &6i6 3 0i + vg\0 2 9 3 9 4 + v$k0 3 94.0 s 

92(0) =^010304 + U^ 5 010405 + vi%e 2 e 3 9 4 + ^45^20405 

g 3 (9) =v^6x9i6i + v[%9 ± 9 3 9 4 + v^ 5 9 2 9 3 9 5 

54(0) =^23010203 + vi&O&ei + ^010405 + W&sfcMfi + B^fcff^ 
5s (0) =^3 W 3 + 010205 + «i 5 3 4 010304 + ugjjWs, 



where the coefficients satisfy v^ k \ = 1 for any triple {j, k, 1} appearing on the right- 
hand-side in the formulas above. The rational map ^ is therefore given by 



9 

P - 



where the sum in the denominator extends over all distinct triples {j,k,l} excluding 
{1,2,4} and {1,3,5}, i.e., the ones corresponding to degenerate cones. 



2.2 Multilinear representation 

The model representation via a rational map ([2]) may be equivalently described in 
terms of a simpler polynomial map ([l]) as follows. Let us substitute 0^ = 9 i s~ 1 ^ d for 
i = 1 , . . . m and define 
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Note that gi : R™ — > R" and Z(0) = Thus we may consider a following more 
convenient version of Q. Find 

6* = argmaxg 1(9) 

n 

subject to Y,h(i)'--9*(d) =^9i{9) = = h V, 9 t > 0. (*) 

CJ" 2 2=1 

Consider a fixed d-tuple of reactions (say, o~\) and in the formulas for gi (i = 1 . . . , n) 
substitute ■ ■ ■ 0«n(d) = 1 ~ So-^o-i ^o-(i) ' ' ' @<r(d)- Note that the resulting algebraic 

statistical map is multilinear i.e, linear in one parameter 9k when all others are fixed. 
For instance, as a function of 9\ we have 

Pi{Q\\') = a$\ + h i = l...,n 

where J2i a i = an d Y^,i bi = l and a,, 6, are given in terms of 8i for Z > 2. 

By Varchenko's theorem (see, [12] chapter 1) the conditional, one dimensional ver- 
sion of problem (*) may be now solved iteratively for each Pi(9i\-), i = l,...,n by 
finding a unique root of the score equations in the regions bounded by the ratios 
—bi/oi. 

Maximization algorithm. Due to the conditional convexity of the one dimensional 
problems the above considerations suggest that the following algorithm for (local) 
maximization of 1(9) should be valid (cf. also [12], Example 1.7, page 11): 

Algorithm 2.1. 

1. Pick initial vectors 9 and 9 Q id € K m . 

2. While \l(9)-l(9 id)\ > e 

• Qold 

• for k=l to m do 

— compute di,bi ( as functions of 9j , j =/= k) 

— identify the bounded interval as determined by Varchenko 's fromula which 
is statistically meaningful (there is only one). 

— use a simple hill-climbing algorithm to find an optimal 9° k vt in that in- 
terval 

— update 9k <— 9° k vt 

3. Recover 9 from 9 by taking 9k = 9k/ @i ■ 

The advantage of the algorithm above is that it reduces a potentially very com- 
plicated multivariate optimization problem in which d and m are large to iteratively 
solving of a simple, univariate one. The disadvantage is that due to its dimension- 
iterative character the algorithm is seen to be slow and for smaller networks perhaps 
less efficient than some off-the-shelf optimization algorithms available in commercial 
software (e.g., some modified hill-climbing methods with random restarts). For that 
reason in our numerical example below we used the standard Matlab optimization 
package rather than Alg. |2.1| 
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In the reminder of the paper we revert to the notation of Section 1 and the original 
problem Q. Based on (*) in this section we may thus extend map g to R> , take 
s(9) = 1 in ([3]) and re-cast the original likelihood maximization problem Q as 

6 = argmaxg^ujogg.^) subject to ^ flq-(i) ■ ■ ■ 6 a [ d ) = 1 and d { > |4j) 
where the g^s are given by ([lj. 

3 Simulated Numerical Example 

In this section we illustrate the ideas discussed above by analyzing a specific numerical 
example in detail. 

If we have d chemical species and data of the form T> — {Ki, i — 1 . . . , k}, then we 
would hope that the statistical algorithm described above should recover the most likely 
d reactions out of a given list of m > d possible reactions, by finding the maximizing 
vector 9 of the corresponding log-likelihood function. In what follows the setup of the 
problem is that of (|4j). To this end, consider the following four-dimensional example: 

M 2A 3 




A l + A 2 A 2 + A 3 

where Ai, i — 0, 1,2,3, denote four chemical species. We shall use the above reaction 
network to simulate "experimentally measured" data and to test the performance of 
our method outlined in Section 2. To this end we shall augment the above network by 
including one or more "incorrect" reactions, and shall check whether our likelihood- 
based algorithm (|4]) is able to identify the original "correct" set of four reactions. 

Data generation. Note that the (deterministic) dynamics of the chemical reaction 
network ([5]) is governed by the linear differential equations of the form 

dAi/dt = jiA i = 0,...,3. (6) 

Thus in our example each data point Ki S ~D (i = 1, . . . , k) was generated by estimating 
the set of parameters (7^) of the true reaction network For each one of the k data 
points the set of parameters (7.;) was drawn independently from a gamma distribution 
G(a, A) with parameters a = 1.5 and A = 1. In order to identify the coordinates 
of the points in T>, the estimated parameters i — 0,1,2,3, were calculated each 
time by fitting the trajectories ^ to the time series data points generated from the 
stochastic process tracing Q (see [S]). The Gillespie algorithm (see [13]) was used 
to generate the 20 equally-spaced values of the trajectory of random process on the 
interval (0, 1) with the fixed initial condition. An example of three random trajectories 
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with independently generated reaction constants values is given in Figure [2] These 
three trajectories would give rise to three independently estimated sets of values (7^) 
and consequently to three data points Ki 6 T>. The fitting was based on the least- 
squares criterion which is statistically justified for estimation purpose of (74) in this 
particular case by an appropriate central limit theorem (cf. e.g., [5] chapter 11). In 
our example, rather than the conditional Algorithm |2.1| we have implemented a more 
widely used local optimization procedure with random restarts as offered by the Matlab 
function lsqcurvef it. 

A resulting single data point Ki = 7^ is the statistical least-squares estimate of 
a realization of four independent gamma-variates G(1.5, 1) which may be viewed as 
coordinates of the true reaction constants vector in the species coordinate system 




Species 
AO-e- 
A1 -e- 



A2 
A3 



0-0- 



.q-O-O-O-0-0-0 



.0 a * = 



\ Q \o A','* A ,Ac6=$r 



O O -0--0- O O O-O 




O A J9-^ A -A-A-A- A-A-A-A-A-A-A- 



!a'''o '' \ ft - - 8.- 0-0-0-0-0-0-0-0-0-0 



0.0 



0.2 



- r— 

0.4 



-0-0 
— I 

0.6 



■0-0-0-0-0-0- 



0.8 



1.0 



time 



Figure 2: An example of generation of the data points Ki € T> for i = 1,2,3 via a two-step 
process of simulation and estimation. Three stochastic trajectories of the reaction network ^ were 
simulated via Gillespie algorithm with propensity (reaction) constants drawn randomly according 
to gamma G(1.5, 1) distribution. The trajectories values at the data collection points are marked at 
20 equally-spaced time-points from to 1. The data from the set of trajectories was used in order 
to estimate the coordinates for each set of coordinates of Ki, i = 1,2,3. The numerical values of 
the reaction rates corresponding to the given trajectories along with their least-squares estimates 
are presented in Table [3j 
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Reaction 


Propensity Constants 


Estimated Values 


Estimators SEs 


A -» A 2 + A 3 
A -> Ai 
4) -> A x + A 2 
4 -> 2A 3 


(0.953, 0.630, 0.065) 
(2.982, 1.869, 0.711) 
(0.328, 0.336, 1.740) 
(1.996, 1.262, 0.853) 


(0.898, 0.668, 0.056) 
(2.711, 1.905, 0.876) 
(0.301, 0.349, 1.874) 
(2.064, 1.212, 0.824) 


(0.091, 0.066, 0.032) 
(0.136, 0.098, 0.060) 
(0.058, 0.048, 0.087) 
(0.118, 0.075, 0.058) 



Table 1: Three sets of reaction kinetic constants corresponding to the trajectories depicted in 
Figure [2] along with their estimated values (obtained via least-squares fitting) and standard errors 
of the estimates. 



[Aq, Ai, A 2 , A 3 ]. This representation of data points is not related to a choice of the 
reactions; note, however, that each data poinlj^jlies inside the open convex cone gener- 
ated by the true reactions ( 5 ) . As shown in pQ the coordinates of Ki in the basis given 
by the reaction vectors in (5 1 are precisely the estimates of the true rate constants. 

The data set V = {Ki,i = 1 . . . , k} used in the simulation described above was 
based on considered k — 50 data points. The first three data points are summarized 
in Table |3 

In order to test our method, let us first add one incorrect reaction, Aq — ► A2, 
and from this point on suppose we have no a priori knowledge of the true chemistry; 
therefore, the five possible reactions are as follows: 




(7) 



Later in this section we also consider the case where we add not just one, but several 
incorrect reactions. 



Calculation of the log-likelihood function. In order to obtain an estimate 9 via 
|4f) one needs to be able to evaluate the map ([TJ, i.e., in addition to the data counts 
vector u e W l in Q one also needs to know the values of the coefficients of the poly- 
nomial map. Whereas the calculation of the exact values is difficult for d > 2, one may 
typically resort to Monte-Carlo approximations (see, e.g., [9]). In our current example, 
for a non-degenerate (i.e., 4-dimensional) cone C, we have computed the approximate 
relative volumes vol(C H Si) /vol(C) using the following Monte Carlo method. For each 

2 Here we assume tacitly that the estimation error is sufficiently small and that the statistical estimation 
procedure is consistent. It turns out this is typically the case in the settings similar to our simulated 
example, but the discussion of the precise conditions under which this is true in real experimental settings 
goes beyond the scope of our present discussion. For our current example a brief inspection of the Table [3] 
indicates a reasonably good agreement between the estimates and the true values of the reaction rates both 
in terms of actual values as well as the corresponding SE's. 
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cone C we generated N — 2000 points inside C with the corresponding conical coor- 
dinates randomly drawn from the four independent gammaG(1.5, 1) random variable 
and then counted the proportion of the total points falling into C n Si i.e., used the 
approximation 



vol(CnSi) 
vol(C) 



(#points mCf]Si)/N i = l...,n. 



With the coefficient values determined as above, the coordinate polynomial maps gi in 
([!]) were easily calculated now by identifying the cones that contained the appropriate 
building block regions Si. 





Figure 3: Geometry of building blocks for 
reaction network ([7]). 



Figure 4: Faces of polyhedron corresponding 
to reaction network ([7]). 



Visualization of the chemical network. The geometry in ^ can be visualized 
in the 3-dimensional subspace W C K n generated by {A 1: A 2 , A 3 }. This follows as all 
the reaction targets are in this subspace, and we can understand the configuration of 
relevant four-dimensional cones by looking at their intersections with W. Each four- 
dimensional cone with vertex Xo intersects W along a tetrahedron. The intersections 
of all these tetrahedra cut out the building blocks corresponding to our example ^ , as 
illustrated in Figures [3] and [4j There are five vertices labeled by numbers corresponding 
to the five target reactions in |7]); they form a six-faced convex polyhedron V. Let C 
be the intersection of line passing through points 1 and 2, denoted (12), with the plane 
(345). Then all building blocks are tetrahedra with a vertex at C and the opposite 
face being one of the six faces of the polyhedron V . For example, the building block 
(C245) is depicted in Figure [3j 

Not surprisingly, the 50 data points generated in our example were distributed 
among the building blocks that compose the tetrahedron (1234) corresponding to the 
true reactions; 32 data points fell inside the building block (C234) and 18 inside (C134). 
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The log-likclihood function was found in this case as 



1(9) = 321og(.706 • 0x020304 + .35 • 9 2 9 3 9 4 9 5 ) + 181og(.294 • 9 1 9 2 9 3 9 i + .339 ■ 0i0 3 4 5 ). 



Maximization of log-likelihood. In order to maximize 1(9) or, equivalently, to 
minimize —1(9), we used the Matlab function fmincon for constrained optimization. 
As in (4J), the constrain is given by the condition s(9) = 9 a (i) 1 ■ ■ 0<r(4) = 1 an( i 
comes from the fact that there are 4 reactions in the true network. This constrain also 
assures that g maps into the probability simplex A„_i, i.e., defines an algebraic variety 
(polynomial map) which corresponds to a valid tstatistical model. 

The optimization is repeated 2 m_1 times (i.e., 16 times for the example for network 
(|7|) with random initial conditions satisfying the constrain. A list of (local) minima 
was created and entries were merged if they were very close. The point 9 that achieved 
the smallest local minimum was reported together with the percentage of time the 
algorithm ended up at that particular point (success rate). The output for example 
([7]) given by the customized Matlab function was 

Minimum of negative log-likelihood: 33.19 
Theta: 

11110. 
Hits: 16 out of 16, 1007. . 

As we may see from these results, in the notation of Figures [3] and [4] the algorithm 
identified the true reactions (targets) 1, 2, 3 and 4 and discarded the incorrect reac- 
tion 5. 

More numerical comparisons. We also ran example ^ with, respectively, two, 
three, four and five incorrect reactions added to the set of four correct ones. The 
true network was always identified and the success rate (percentage of correct hits for 
various random initial guess) was high. The results of these additional experiments are 
summarized in Table 2. 



# reactions 


# cones 


# non-degenerate 


# building 


running 


success 






cones 


blocks 


time 


rate 


m=5 


5 


5 


6 


7 sec 


100% 


m=6 


15 


11 


15 


29 sec 


94% 


m=7 


35 


30 


133 


6 min 32 sec 


94% 


m=8 


70 


64 


871 


lh 12 min 


98% 


m=9 


126 


115 


2397 


8h 55 min 


96% 



Table 2: Summary of numerical results for k — 50 data points, using N = 2000 rays in the Monte 
Carlo relative volume computation, 2 m ~ 1 optimizations with random initial guess. The code was 
tested on a 2.8 Ghz Intel Core2Duo iMac machine. 
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4 Summary and Discussion 



We have proposed herein a statistical method for inferring a biochemical reaction net- 
work given several sets of data that originate from "noisy" versions of the reaction rate 
equations associated with the network. As illustrated in some earlier work of some of 
the authors pQ, in the usual deterministic sense such networks are in general uniden- 
tifiable, i.e., different chemical reaction networks may give rise to exactly the same 
reaction rate equations. In practice, the matters are further complicated since the co- 
efficients of the reaction rate equations are estimated from available experimental data, 
and hence are subject to measurement error and, moreover, their actual values may 
differ at different experimental conditions, i.e. at different data points. The statistical 
approach described here is largely unaffected by these problems, as it only relies on 
the geometry of the network relative to the data distribution, in order to identify the 
sets of most likely reactions. Hence, the method takes advantage of the algebraic and 
geometric representation of the network rather than merely the observed experimen- 
tal values of the network species, as is commonly the case in network inference models 
based on graphical methods, like e.g. Baycsian or probabilistic boolean networks. Still, 
in order to use the proposed multinomial parametrization of a biochemical network, 
the method does require a valid way of mapping the experimentally estimated rate 
coefficients into the networks' appropriate convex regions, and with very large mea- 
surement errors is likely to perform poorly. On the other hand, precisely because of 
the need for the experimental data mapping, the method has a very attractive feature 
of being able to potentially combine variety of different data sets obtained by various 
methods into one set of experimental points placed in a convex hull of the network 
building blocks. These universality properties of the method require further studies 
and possibly a development of additional statistical methodology beyond the scope of 
our present work. In the current paper our main goal was to present a proof-of-concept 
example based on simulated data, with a purposefully straightforward but non-trivial 
model discrimination problem. For the example provided in this paper the method was 
seen to perform very well, with almost perfect discrimination against incorrect models 
even as the complexity of the model selection problem increased. 

Nonetheless, further studies and developments are needed to assess how well the 
method may perform on more challenging and realistic data sets. In particular, one of 
the aspects of the methodology which was not pursued here, and which could improve 
its computational scalability, is the utilization of techniques from computational algebra 
in order to increase the efficiency and further automate the proposed maximization 
algorithm. 
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